home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Support / Convolution.m next >
Encoding:
Text File  |  1992-08-18  |  24.9 KB  |  713 lines
  1. (*  :Title:    Piecewise Convolution Rule Base  *)
  2.  
  3. (*  :Authors:    Kevin West and Brian Evans  *)
  4.  
  5. (*
  6.     :Summary:    To take the continuous-time and discrete-time
  7.         convolution of two piecewise-defined functions and return
  8.         another piecewise-defined function.
  9.  *)
  10.  
  11. (*  :Context:    SignalProcessing`Support`Convolution`  *)
  12.  
  13. (*  :PackageVersion:  2.7    *)
  14.  
  15. (*
  16.     :Copyright:    Copyright 1989, 1990 by
  17.         the Digital Signal Processing Group
  18.         at the Georgia Institute of Technology
  19.  
  20.     Permission to use, copy, modify, and distribute this software
  21.     and its documentation for any purpose and without fee is
  22.     hereby granted, provided that the above copyright notice
  23.     appear in all copies and that both that copyright notice and
  24.     this permission notice appear in supporting documentation,
  25.     and that the name of Georgia Tech or Georgia Institute of
  26.     Technology not be used in advertising or publicity pertaining
  27.     to distribution of the software without specific, written prior
  28.     permission.  Georgia Tech makes no representations about the
  29.     suitability of this software for any purpose.  It is provided
  30.     "as is" without express or implied warranty.
  31.  *)
  32.  
  33. (*  :History:    *)
  34.  
  35. (*  :Keywords:    *)
  36.  
  37. (*  :Source:    {Discrete-Time Signal Processing} by Oppenheim & Schafer  *)
  38.  
  39. (*  :Warning:    *)
  40.  
  41. (*  :Mathematica Version:  1.2 or 2.0  *)
  42.  
  43. (*  :Limitation:  *)
  44.  
  45. (*
  46.     :Discussion:  This package originally performed piecewise continuous
  47.           convolution and therefore was part of the analog signal
  48.           processing packages (spp).  We've added discrete convolution,
  49.           so it is now in the Support packages.  It's old context
  50.           was "SignalProcessing`Analog`Convolution`"---  if you
  51.           ask Mathematica to load Convolution from the analog spp,
  52.           it will still load.
  53.  *)
  54.  
  55. (*
  56.     :Functions:    AutoCorrelation
  57.         ConvertFromList
  58.         ConvertToList
  59.         CTPiecewiseConvolution
  60.         DTPiecewiseConvolution
  61.         IntervalQ
  62.         PiecewiseConvolution
  63.         PlotList
  64.         SimplifyList
  65.         ValidIntervalQ
  66.  *)
  67.  
  68.  
  69.  
  70. (*  B E G I N     P A C K A G E  *)
  71.  
  72. BeginPackage [    "SignalProcessing`Support`Convolution`",
  73.         "SignalProcessing`Support`SigProc`",
  74.         "SignalProcessing`Support`SupCode`",
  75.         If [ TrueQ[$VersionNumber >= 2.0],
  76.              "Algebra`SymbolicSum`",
  77.              "Algebra`GosperSum`" ] ]
  78.  
  79.  
  80. If [ TrueQ[ $VersionNumber >= 2.0 ],
  81.      Off[ General::spell ];
  82.      Off[ General::spell1 ] ]
  83.  
  84.             
  85. (*  U S A G E     I N F O R M A T I O N  *)
  86.  
  87. $ConvolutionDomain::usage = 
  88.     "$ConvolutionDomain is the default domain for all of the routines \
  89.     in the Convolution package (except for DTPiecewiseConvolution and \
  90.     CTPiecewiseConvolution which imply a domain). \
  91.     Use SetConvolutionDomain to reset its value."
  92.  
  93. Area::usage =
  94.     "Area[ar] is an object representing the area under a Dirac delta \
  95.     function. \
  96.     A Dirac delta function is usually written as C Delta[t - t0], \
  97.     where C is the area under the delta function and t0 is the location \
  98.     of the delta function. \
  99.     In F-interval form, C Delta[t - t0] becomes { Area[C], t0, t0 }."
  100.  
  101. AutoCorrelation::usage =
  102.     "AutoCorrelation[e, v] computes the autocorrelation of expression e \
  103.     (in the form of a list of F-intervals, an F-interval, or a function) \
  104.     with respect to v. \
  105.     The variable v is treated as either discrete and continuous \
  106.     according to the value of the Domain option. \
  107.     The default domain is the value of $ConvolutionDomain."
  108.  
  109. ConvertFromList::usage =
  110.     "ConvertFromList[l, t] converts a list of F-intervals in the form \
  111.     {f, t1, t2} to a signal processing expression. \
  112.     The new expression will contain either Delta, CStep, and CPulse \
  113.     functions (continuous-time) or Impulse, Step, and Pulse functions \
  114.     (discrete-time) according to the value of the Domain option. \
  115.     The default domain is the value of $ConvolutionDomain."
  116.  
  117. ConvertToList::usage =
  118.     "ConvertToList[e, t] converts the expression e to a \
  119.     list of F-intervals in the form {fun, left, right}. \
  120.     The domain of the new expression is established according to the \
  121.     value of the Domain option. \
  122.     The default domain is the value of $ConvolutionDomain."
  123.  
  124. CTPiecewiseConvolution::usage =
  125.     "CTPiecewiseConvolution[f, g, t] carries out the one-dimensional,
  126.     continuous-time piecewise convolution of f and g with respect to t. \
  127.     See PiecewiseConvolution for valid representation of f and g. \
  128.     Using F-interval notation {function, left_endpoint, right_endpoint}, \
  129.     a Dirac delta is meant when function = Area[<area>] \
  130.     (where <area> equals the area of under the impulse function) and \
  131.     left_endpoint = right_endpoint (which is the location of the impulse)."
  132.  
  133. DTPiecewiseConvolution::usage = 
  134.     "DTPiecewiseConvolution[f, g, n] carries out the one-dimensional,
  135.     discrete-time piecewise convolution of f and g with respect to n. \
  136.     See PiecewiseConvolution for the valid representations of f and g. \
  137.     Using F-interval notation {function, left_endpoint, right_endpoint}, \
  138.     an impulse (Kronecker delta) function is meant when function = <value> \
  139.     (where <value> equals the strength of the impulse function) and \
  140.     left_endpoint = right_endpoint (which is the location of the impulse)."
  141.  
  142. IntervalQ::usage =
  143.     "IntervalQ[x] returns True if the argument x is an interval, and \
  144.     gives False otherwise. \
  145.     An interval has the form {f, t1, t2} such that none of f, t1, or t2 \
  146.     are lists. \
  147.     See also ValidIntervalQ."
  148.  
  149. PiecewiseConvolution::usage =
  150.     "PiecewiseConvolution[f, g, v] convolves the piecewise functions \
  151.     f and g with respect to the variable v. \
  152.     In this context, a function is represented in a piecewise fashion: \
  153.     (1) as F-intervals of the form {fun, left, right}, \
  154.     (2) as a list of F-intervals, or (3) as an expression. \
  155.     An F-interval has the form {function, left_endpoint, right_endpoint}. \
  156.     The F-interval notation represents a finite-extent function or \
  157.     sequence when the endpoints do not equal infinity. \
  158.     The result of convolution is always returned as a list of F-intervals. \
  159.     The variable v is treated as either discrete and continuous \
  160.     according to the value of the Domain option. \
  161.     The default domain is the value of $ConvolutionDomain, \
  162.     which can be reset by SetConvolutionDomain. \
  163.     See also CTPiecewiseConvolution and DTPiecewiseConvolution."
  164.  
  165. PlotList::usage =
  166.     "PlotList[l, {v, vstart, vend}] plots the piecewise function l vs. v. \
  167.     The variable v is treated as either discrete and continuous \
  168.     according to the value of the Domain option. \
  169.     The default domain is the value of $ConvolutionDomain."
  170.  
  171. SetConvolutionDomain::usage =
  172.     "SetConvolutionDomain[domain] establishes the default convolution \
  173.     domain as either Continuous or Discrete. \
  174.     It returns domain on success, and Null on error."
  175.  
  176. SimplifyList::usage =
  177.     "SimplifyList[interval_list, v] takes a list of F-intervals \
  178.     each having the form {f, t1, t2} and returns a simplified sorted \
  179.     list of F-intervals. \
  180.     The interval_list is treated as either a continuous-time or a \
  181.     discrete-time function according the value of the Domain option. \
  182.     The variable v is treated as either discrete and continuous \
  183.     according to the value of the Domain option. \
  184.     The default domain is the value of $ConvolutionDomain."
  185.  
  186. ValidIntervalQ::usage =
  187.     "ValidIntervalQ[i] returns True if i is a valid F-interval, and \
  188.     gives False otherwise. \
  189.     An F-interval is an interval of the form {f, leftendp, rightendp} \
  190.     where f is a function and leftendp < rightendp unless f is an \
  191.     impulse function in which case leftendp must equal rightendp. \
  192.     In a valid F-interval, leftendp and rightenp must be numbers. \
  193.     See also IntervalQ."
  194.  
  195. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  196.  
  197.  
  198.  
  199. Begin[ "`Private`" ]
  200.  
  201.  
  202. (*  S Y S T E M     S E T T I N G S  *)
  203.  
  204. $ConvolutionDomain = Continuous
  205.  
  206. Protect[ $ConvolutionDomain ]
  207.  
  208.  
  209. (*  M E S S A G E S  *)
  210.  
  211. AutoCorrelation::domain =
  212.     "The option Domain must be equal to Discrete or Continuous."
  213. ConvertFromList::domain =
  214.     "The option Domain must be equal to Discrete or Continuous."
  215. ConvertToList::domain =
  216.     "The option Domain must be equal to Discrete or Continuous."
  217. CTPiecewiseConvolution::argct =
  218.     "CTPiecewiseConvolution requires exactly three arguments."
  219. CTPiecewiseConvolution::badint =
  220.     "The interval `` is not valid and has been deleted."
  221. DTPiecewiseConvolution::argct =
  222.     "DTPiecewiseConvolution requires exactly three arguments."
  223. DTPlotList::badargs =
  224.     "DTPlotList requires the format DTPlotList[list,{n,n1,n2}]"
  225. PiecewiseConvolution::argct =
  226.     "PiecewiseConvolution requires at least three arguments."
  227. PiecewiseConvolution::domain =
  228.     "The option Domain must be equal to Discrete or Continuous."
  229. PlotList::badargs =
  230.     "PlotList requires the format PlotList[list,{t,t1,t2}]"
  231. PlotList::domain =
  232.     "The option Domain must be equal to Discrete or Continuous."
  233. SimplifyList::domain =
  234.     "The option Domain must be equal to Discrete or Continuous."
  235. SetConvolutionDomain::domain =
  236.     "The option Domain must be equal to Discrete or Continuous."
  237. ValidIntervalQ::domain =
  238.     "The option Domain must be equal to Discrete or Continuous."
  239.  
  240. (*  E N D     M E S S A G E S  *)
  241.  
  242.  
  243. (*  S U P P O R T I N G     R O U T I N E S  *)
  244.  
  245.  
  246. (*  AreaOfImpulse -- return the impulse area at a particular point *)
  247. AreaOfImpulse[list_, a_] :=
  248.    Block[ {areaofimp = 0, i},
  249.           Do[ If [ Head[list[[i]][[1]]] == Area,
  250.                    areaofimp += list[[i]][[1]][[1]] * ImpulseHere[list[[i]], a],
  251.                    Null ],
  252.               {i, Length[list]} ];
  253.           areaofimp ]
  254.  
  255. (*  CallFunction  *)
  256. SetAttributes[CallFunction, {HoldFirst}]
  257.  
  258. CallFunction[ message_, discretefun_, continuousfun_, oplist_List, args___ ] :=
  259.     Block [    {domain},
  260.         domain = Replace[Domain, oplist];
  261.         Which [    SameQ[domain, Discrete],
  262.               Apply[discretefun, {args}],
  263.             SameQ[domain, Continuous],
  264.               Apply[continuousfun, {args}],
  265.             True,
  266.               message ] ]
  267.  
  268.  
  269. (*  ConvertToList routines  *)
  270.  
  271. Options[ConvertToList] := { Domain :> $ConvolutionDomain }
  272.  
  273. ConvertToList[f_, v_, options___] :=
  274.     CallFunction[ Message[ConvertToList::domain],
  275.               DTConvertToList, CTConvertToList,
  276.               ToList[options] ~Join~ Options[ConvertToList],
  277.               f, v ]
  278.  
  279. CTConvertToList[f_, t_] := CTSimplifyList[ UnsimpCTL[f, t], t ]
  280. DTConvertToList[f_, n_] := DTSimplifyList[ UnsimpCTL[f, n], n ]
  281.  
  282. UnsimpCTL[f_?IntervalQ, t_] := {f}
  283. UnsimpCTL[f_List, t_] := f
  284. UnsimpCTL[f_ + g_, t_] := Join[UnsimpCTL[f,t], UnsimpCTL[g,t]]
  285. UnsimpCTL[0, t_] := {}
  286. UnsimpCTL[Delta[t_ + a_.] ar_., t_] :=
  287.     {{Area[ Limit[ar, t -> -a] ], -a, -a}} /; FreeQ[a,t]
  288. UnsimpCTL[Delta[-t_ + a_.] ar_., t_] :=
  289.     {{Area[ Limit[ar, t -> a] ], a, a}} /; FreeQ[a,t]
  290. UnsimpCTL[Impulse[t_ + a_.] ar_., t_] :=
  291.     {{Limit[ ar, t -> -a], -a, -a}} /; FreeQ[a,t]
  292. UnsimpCTL[Impulse[-t_ + a_.] ar_., t_] :=
  293.     {{Limit[ar, t -> a], a, a}} /; FreeQ[a,t]
  294. UnsimpCTL[CStep[t_ + a_.] f_., t_] := {{f,-a,Infinity}} /; FreeQ[a, t]
  295. UnsimpCTL[CStep[-t_ + a_.] f_., t_] := {{f,-Infinity,a}} /; FreeQ[a, t]
  296. UnsimpCTL[Step[t_ + a_.] f_., t_] := {{f,-a,Infinity}} /; FreeQ[a, t]
  297. UnsimpCTL[Step[-t_ + a_.] f_., t_] := {{f,-Infinity,a}} /; FreeQ[a, t]
  298. UnsimpCTL[CPulse[l_, t_ + a_.] f_.,t_] := {{f, -a, l - a}} /; FreeQ[{l,a}, t]
  299. UnsimpCTL[CPulse[l_, -t_ + a_.] f_.,t_] := {{f, a - l, a}} /; FreeQ[{l,a}, t]
  300. UnsimpCTL[Pulse[l_, t_ + a_.] f_.,t_] := {{f, -a, l - 1 - a}} /; FreeQ[{l,a}, t]
  301. UnsimpCTL[Pulse[l_, -t_ + a_.] f_.,t_] := {{f, a - l + 1, a}} /; FreeQ[{l,a}, t]
  302. UnsimpCTL[f_ term_Plus, t_] := UnsimpCTL[ Distribute[f term], t ] /;
  303.     ! FreeQ[term, t] && ! SameQ[f term, Distribute[f term]]
  304. UnsimpCTL[f_, t_] := {{f,-Infinity,Infinity}}
  305.  
  306. (*  ConvertFromList routines  *)
  307.  
  308. Options[ConvertFromList] := { Domain :> $ConvolutionDomain }
  309.  
  310. ConvertFromList[list_, v_, options___] :=
  311.     CallFunction[ Message[ConvertFromList::domain],
  312.               DTConvertFromList, CTConvertFromList,
  313.               ToList[options] ~Join~ Options[ConvertFromList],
  314.               list, v ]
  315.  
  316. CTConvertFromList[list_, t_] :=
  317.     Apply[Plus, Map[CTConvertFromInterval[#1, t]&, list]]
  318. DTConvertFromList[list_, t_] :=
  319.         Apply[Plus, Map[DTConvertFromInterval[#1, t]&, list]]
  320.  
  321. CTConvertFromInterval[{Area[ar_], a_, a_}, t_] := ar Delta[t - a]
  322. CTConvertFromInterval[{f_, -Infinity, Infinity}, t_] := f
  323. CTConvertFromInterval[{f_, a_, Infinity}, t_] := f CStep[t - a]
  324. CTConvertFromInterval[{f_, -Infinity, a_}, t_] := f CStep[a - t]
  325. CTConvertFromInterval[{f_, a_, b_}, t_] := f CPulse[b - a, t - a]
  326.  
  327. DTConvertFromInterval[{Area[ar_], a_, a_}, t_] := ar Impulse[t-a]  
  328. DTConvertFromInterval[{ar_, a_, a_}, t_] := ar Impulse[t-a]
  329. DTConvertFromInterval[{f_, -Infinity, Infinity}, t_] := f
  330. DTConvertFromInterval[{f_, a_, Infinity}, t_] := f Step[t - a]
  331. DTConvertFromInterval[{f_, -Infinity, a_}, t_] := f Step[a - t]
  332. DTConvertFromInterval[{f_, a_, b_}, t_] := f Pulse[b - a + 1, t - a]
  333.  
  334.  
  335. (*  ConvolveTwoIntervals  *)
  336. ConvolveTwoIntervals[pair_List, t_, time_] :=
  337.     ConvolveTwoIntervals[ pair[[1]], pair[[2]], t, time]
  338.  
  339. ConvolveTwoIntervals[{Area[ar1_], a_, a_}, {Area[ar2_], b_, b_}, t_, time_] :=
  340.     {{Area[ar1 ar2], a+b, a+b}} 
  341. ConvolveTwoIntervals[{Area[ar_], a_, a_}, {f_, b_, c_}, t_, time_] :=
  342.     {{ar (f /. t -> t-a), b+a, c+a}}
  343. ConvolveTwoIntervals[{f_, a_, b_}, {Area[ar_], c_, c_}, t_, time_] :=
  344.     {{ar (f /. t -> t-c), a+c, b+c}}
  345. ConvolveTwoIntervals[{f_, -Infinity, t2_}, {g_, t3_, Infinity}, t_, time_] :=
  346.     ConvolveTwoIntervals[{g, t3, Infinity}, {f, -Infinity, t2}, t, time] /;
  347.     t3 > -Infinity
  348.       (* Jan-92 JMc: make sure left-sided signal is g(t) *)
  349. ConvolveTwoIntervals[{f_, t1_, t2_}, {g_, t3_, t4_}, t_, time_] :=
  350.     ConvolveTwoIntervals[{g, t3, t4}, {f, t1, t2}, t, time] /;
  351.     (t4 - t3) > (t2 - t1)
  352.  
  353. ConvolveTwoIntervals[{f_, t1_, t2_}, {g_, t3_, t4_}, t_, Continuous] := 
  354.     Block [    {convlist, f2, g2, T, ta, tb, tc, td},
  355.         f2 = f /. t -> T;
  356.         g2 = g /. t -> t - T;
  357.         ta = If [ t1 > -Infinity, t1+t4,  t1 ];
  358.         tb = If [ t2 <  Infinity, t2+t3,  t2 ];
  359.         tc = If [ t3 > -Infinity, t -t3, -t3 ];
  360.         td = If [ t4 <  Infinity, t -t4, -t4 ];
  361.         convlist = {};
  362.         If [ ta > -Infinity,
  363.              AppendTo[convlist,
  364.                       {Integrate[f2 g2, {T,t1,tc}], t1+t3, ta}] ];
  365.         If [ ta < Infinity && tb > -Infinity,
  366.              AppendTo[convlist,
  367.                       {Integrate[f2 g2, {T,td,tc}], ta, tb}] ];
  368.         If [ tb < Infinity,
  369.              AppendTo[convlist,
  370.                       {Integrate[f2 g2, {T,td,t2}], tb, t2+t4}] ];
  371.         convlist ]
  372.  
  373. ConvolveTwoIntervals[{f_, t1_, t2_}, {g_, t3_, t4_}, t_, Discrete] := 
  374.     Block [    {convlist, f2, g2, sumhead, summation, tt1, tt2, T, ta, tb, tc, td},
  375.         f2 = f /. t -> T;
  376.         g2 = g /. t -> t - T;
  377.         ta = If [ t1 > -Infinity, t1+t4,  t1 ];
  378.         tb = If [ t2 <  Infinity, t2+t3,  t2 ];
  379.         tc = If [ t3 > -Infinity, t -t3, -t3 ];
  380.         td = If [ t4 <  Infinity, t -t4, -t4 ];
  381.         sumhead = If [ TrueQ[$VersionNumber >= 2.0],
  382.                    SymbolicSum,
  383.                    GosperSum ];
  384.         summation = Cancel[sumhead[f2 g2, {T, tt1, tt2}]];
  385.         convlist = {};
  386.         If [ ta > -Infinity,
  387.              AppendTo[convlist,
  388.                       {summation /. tt1 -> t1 /. tt2 -> tc, t1+t3, ta}] ];
  389.         If [ ta < Infinity && tb > -Infinity,
  390.              AppendTo[convlist,
  391.                       {summation /. tt1 -> td /. tt2 -> tc, ta+1, tb}] ];
  392.         If [ tb < Infinity,
  393.              AppendTo[convlist,
  394.                       {summation /. tt1 -> td /. tt2 -> t2, tb+1, t2+t4+1}] ];
  395.         convlist ]
  396.  
  397. (*  ConvTwoLists  *)
  398. PairIntervalAndList[int_, list_] := Map[ {int, #}&, list ]
  399.  
  400. CTConvTwoLists[a_, b_, v_] := ConvTwoLists[a, b, v, Continuous]
  401. DTConvTwoLists[a_, b_, v_] := ConvTwoLists[a, b, v, Discrete]
  402.  
  403. ConvTwoLists[{}, b_, v_, time_] := {}
  404. ConvTwoLists[a_, {}, v_, time_] := {}
  405. ConvTwoLists[a_, b_, v_, time_] :=
  406.     Block [    {convolution, groupedinto3s, nestedlist, pairedlists},
  407.         nestedlist = Flatten[ Map[ PairIntervalAndList[#, b]&, a ] ];
  408.         pairedlists = Partition[ Partition[nestedlist, 3], 2 ];
  409.         convolution = Map[ ConvolveTwoIntervals[#, v, time]&,
  410.                    pairedlists ];
  411.         groupedinto3s = Partition[ Flatten[convolution], 3 ];
  412.         SimplifyList[ RemoveBadInterval[groupedinto3s, time], v,
  413.                   Domain -> time ] ]
  414.  
  415. (*  ImpulseHere -- return 1 if the points coincide *)
  416. ImpulseHere[{Area[ar_], a_, a_}, a_] := 1
  417. ImpulseHere[___] := 0
  418.  
  419. (*  InIntervalQ; InInterval returns a 1 if the intervals cross  *)    
  420. InIntervalQ[{f_, left_, right_}, a_, b_] :=
  421.     TrueQ[ N[MyLessEqual[left, a, b, right]] ]
  422. InInterval[args__] := If [ InIntervalQ[args], 1, 0 ]
  423.  
  424. (*  IntervalQ  *)
  425. IntervalQ[{f_, a_, b_}] := ! ( ListQ[f] || ListQ[a] || ListQ[b] )
  426. IntervalQ[___] := False
  427.  
  428. (*  ListOfEndpoints  *)
  429. getsecthird[ {a_, b_, c_} ] := ToCollection[b, c]
  430.  
  431. ListOfEndpoints[list_List] := Sort[ Union[ Map[getsecthird, list] ], MyLess ]
  432.  
  433.  
  434. (*  MyLess  *)
  435. MyLess[a_, a_] := False    
  436. MyLess[-Infinity, Infinity] := True
  437. MyLess[Infinity, -Infinity] := False
  438. MyLess[Infinity, a_] := False
  439. MyLess[a_, Infinity] := True
  440. MyLess[a_, -Infinity] := False
  441. MyLess[-Infinity, a_] := True
  442. MyLess[a_, b_] := N[a] < N[b]
  443. MyLess[a_, b_, rest__] := MyLess[a, b] && MyLess[b, rest]
  444.  
  445. (*  MyLessEqual  *)
  446. MyLessEqual[a_, a_] := True
  447. MyLessEqual[a_, b_] := MyLess[a, b]
  448. MyLessEqual[a_, b_, rest__] := MyLessEqual[a, b] && MyLessEqual[b, rest]
  449.  
  450. (*  RemoveBadInterval  *)
  451. RemoveBadInterval[oldlist_, Discrete] := Select[oldlist, DTValidIntervalQ]
  452. RemoveBadInterval[oldlist_, Continuous] := Select[oldlist, CTValidIntervalQ]
  453.  
  454. (*  RemoveDomain  *)
  455. keepq[x_] := SameQ[Domain, Replace[Domain, x]]
  456. RemoveDomain[oplist_] := Select[oplist, keepq]
  457.  
  458. (*  ValueOfInterval -- return the value of a particular interval  *)
  459. ValueOfInterval[list_, a_, b_] :=
  460.     Apply[Plus, Map[First, Select[list, InIntervalQ[#1, a, b]&]]]
  461.  
  462. (*  ValidIntervalQ  *)
  463. Options[ValidIntervalQ] := { Domain :> $ConvolutionDomain }
  464.  
  465. ValidIntervalQ[list_, options___] := 
  466.     CallFunction[ Message[ValidIntervalQ::domain],
  467.               DTValidIntervalQ, CTValidIntervalQ,
  468.               ToList[options] ~Join~ Options[ValidIntervalQ],
  469.               list ]
  470.  
  471. CTValidIntervalQ[i_] := False /; Not[IntervalQ[i]]
  472. CTValidIntervalQ[{Area[ar_], a_, b_}] := TrueQ[ (a == b) && (ar != 0) ]
  473. CTValidIntervalQ[{f_, a_, b_}] := TrueQ[ N[(b > a) && ! SameQ[f, 0]] ]
  474.  
  475. DTValidIntervalQ[i_] := False /; Not[IntervalQ[i]]
  476. DTValidIntervalQ[{Area[ar_], a_, b_}] :=
  477.     SameQ[a, b] && (! InfinityQ[b]) && (! SameQ[ar, 0])
  478. DTValidIntervalQ[{f_, a_, a_}] := (! InfinityQ[a]) && (! SameQ[f, 0])
  479. DTValidIntervalQ[{f_, a_, b_}] := TrueQ[ N[(b >= a) && (! SameQ[f, 0])] ]
  480.  
  481.  
  482. (*  C O N V O L U T I O N     R O U T I N E S  *)
  483.  
  484.  
  485. (*  AutoCorrelation  *)
  486.  
  487. Options[AutoCorrelation] := { Domain :> $ConvolutionDomain }
  488.  
  489. AutoCorrelation[f_, v_, options___] :=
  490.     CallFunction[ Message[AutoCorrelation::domain],
  491.               DTPiecewiseConvolution, CTPiecewiseConvolution,
  492.               ToList[options] ~Join~ Options[AutoCorrelation],
  493.               f, FlipInTime[f, v], v ]
  494.  
  495. FlipInTime[int_, v_] := FlipInterval[int, v] /; IntervalQ[int]
  496. FlipInTime[list_List, v_] := Reverse[ Map[FlipInterval[#1, v]&, list] ]
  497. FlipInTime[f_, v_] := f /. v -> -v
  498.  
  499. FlipInterval[{f_, v1_, v2_}, v_] := {(f /. v -> -v), -v2, -v1}
  500.  
  501. (*  PiecewiseConvolution  *)
  502. Options[PiecewiseConvolution] := { Domain :> $ConvolutionDomain }
  503.  
  504. PiecewiseConvolution[f_, g_, v_, options___] :=
  505.     CallFunction[ Message[PiecewiseConvolution::domain],
  506.               DTPiecewiseConvolution, CTPiecewiseConvolution,
  507.               ToList[options] ~Join~ Options[PiecewiseConvolution],
  508.               f, g, v ]
  509.  
  510. PiecewiseConvolution[___] := Message[PiecewiseConvolution::argct]
  511.  
  512. Unprotect[CConvolve, Convolve]
  513.  
  514. (*  CTPiecewiseConvolution  *)
  515. CConvolve/: TheFunction[ CConvolve[t_Symbol][x_, y_] ] :=
  516.     CTPiecewiseConvolution[x, y, t]
  517. CConvolve/: TheFunction[ CConvolve[t_Symbol][x_, y_, rest__] ] :=
  518.     TheFunction[ CConvolve[t][ CTPiecewiseConvolution[x, y, t], rest ] ]
  519.  
  520. CTPiecewiseConvolution[f_, g_, t_] :=
  521.     CTConvTwoLists[ CTConvertToList[f, t], CTConvertToList[g, t], t ]
  522. CTPiecewiseConvolution[___] := Message[CTPiecewiseConvolution::argct]
  523.  
  524. (*  DTPiecewiseConvolution  *)
  525. Convolve/: TheFunction[ Convolve[n_Symbol][x_, y_] ] :=
  526.     DTPiecewiseConvolution[x, y, n]
  527. Convolve/: TheFunction[ Convolve[n_Symbol][x_, y_, rest__] ] :=
  528.     TheFunction[ Convolve[n][ DTPiecewiseConvolution[x, y, n], rest ] ]
  529.  
  530. DTPiecewiseConvolution[f_, g_, n_] :=
  531.     DTConvTwoLists[ DTConvertToList[f, n], DTConvertToList[g, n], n ]
  532. DTPiecewiseConvolution[___] := Message[DTPiecewiseConvolution::argct]
  533.  
  534. Protect[CConvolve, Convolve]
  535.  
  536. (*  P L O T T I N G     R O U T I N E S  *)
  537.  
  538. (*  PlotList  *)
  539. Options[PlotList] := { Domain :> $ConvolutionDomain }
  540.  
  541. PlotList[list_, {v_, v1_, v2_}, options___] :=
  542.     Block [    {oplist},
  543.         oplist = ToList[options] ~Join~ Options[PlotList];
  544.         CallFunction[ Message[PlotList::domain],
  545.                   DTPlotList, CTPlotList, oplist,
  546.                   list, {v, v1, v2}, RemoveDomain[oplist] ] ]
  547.  
  548. (* CTPlotList  *)
  549. CTPlotList[int_, {t_, t1_, t2_}, options___] :=
  550.         PlotList[{int}, {t, t1, t2}, options] /; IntervalQ[int]
  551. CTPlotList[list_List, {t_, t1_, t2_}, options___] :=
  552.     SignalPlot[CTConvertFromList[list, t], {t, t1, t2}, options]
  553. CTPlotList[fun_, {t_, t1_, t2_}, options___ ] :=
  554.     SignalPlot[fun, {t, t1, t2}, options]
  555.  
  556. (* DTPlotList *)
  557. DTPlotList[int_, {n_, n1_, n2_}, options___] :=
  558.     DTPlotList[{int}, {n, n1, n2}, options] /; IntervalQ[int]
  559. DTPlotList[list_List, {n_, n1_, n2_}, options___] :=
  560.     SequencePlot[DTConvertFromList[list, n], {n, n1, n2}, options]
  561. DTPlotList[fun_, {n_, n1_, n2_}, options___] :=
  562.     SequencePlot[fun, {n, n1, n2}, options]
  563.  
  564.  
  565. (*  S I M P L I F I C A T I O N     R O U T I N E S  *)
  566.  
  567. (*  SimplifyList  *)
  568. Options[SimplifyList] := { Domain :> $ConvolutionDomain }
  569.  
  570. SimplifyList[list_, v_, options___] :=
  571.     CallFunction[ Message[SimplifyList::domain],
  572.               DTSimplifyList, CTSimplifyList,
  573.               ToList[options] ~Join~ Options[SimplifyList],
  574.               list, v ]
  575.  
  576. (*  CTSimplifyList  *)
  577. CTSimplifyList[{}, t_] := {}
  578. CTSimplifyList[list_, t_] :=
  579.     Block[ { curpoint, k, nextpoint, points, simplelist},
  580.         points = ListOfEndpoints[list];        
  581.         nextpoint = points[[1]];
  582.         simplelist = {{ Area[AreaOfImpulse[list, nextpoint]],
  583.                 nextpoint, nextpoint }};
  584.         Do [ curpoint = nextpoint;
  585.              nextpoint = points[[k]];
  586.              AppendTo[ simplelist,
  587.                    { ValueOfInterval[list, curpoint, nextpoint],
  588.                  curpoint, nextpoint } ];
  589.              AppendTo[ simplelist,
  590.                    { Area[AreaOfImpulse[list, nextpoint]],
  591.                                 nextpoint, nextpoint } ],
  592.              {k, 2, Length[points]} ];
  593.         RemoveBadInterval[SPSimplify[simplelist, Variables -> {t}],
  594.                   Continuous] ]
  595.  
  596. (*  DTSimplifyList  *)
  597. OpenList[{Area[A_], a_, a_}] := {A, a, a+1}
  598. OpenList[{Area[A_], a_, b_}] := {} /; ! SameQ[a, b]
  599. OpenList[{f_, a_, b_}] := {f, a, b+1}
  600. CloseList[{f_, a_, b_}] := {f, a, b-1}
  601.  
  602. DTSimplifyList[{}, n_] := {}
  603. DTSimplifyList[list_, n_] :=
  604.     Block [    {curpoint, ival, k, lastpoint, nextpoint, numpoints,
  605.          openlist, points, simplelist = {}, simplified},
  606.  
  607.         openlist = Map[OpenList, list];
  608.         points = ListOfEndpoints[openlist];
  609.         numpoints = Length[points];
  610.         nextpoint = points[[1]];
  611.         Do [ curpoint = nextpoint;
  612.              nextpoint = points[[k]];
  613.              ival = ValueOfInterval[openlist, curpoint, nextpoint];
  614.              AppendTo[ simplelist, {ival, curpoint, nextpoint} ],
  615.              {k, 2, numpoints} ];
  616.         simplified = SPSimplify[Map[CloseList, simplelist],
  617.                     Variables -> {n}];
  618.         simplelist = RemoveBadInterval[simplified, Discrete];
  619.  
  620.         openlist = {};
  621.         While [    Length[simplelist] > 0,
  622.             curpoint = First[ simplelist ];
  623.             simplelist = Rest[ simplelist ];
  624.  
  625.             If [ ! SameQ[ curpoint[[2]], curpoint[[3]] ],
  626.                  AppendTo[ openlist, curpoint ];
  627.                  curpoint = {},
  628.  
  629.                  If [ Length[openlist] > 0,
  630.                   lastpoint = Last[openlist];
  631.                   If [ SameQ[curpoint[[2]], 1+lastpoint[[3]]] &&
  632.                        SameQ[curpoint[[1]] /. n->curpoint[[2]],
  633.                            lastpoint[[1]] /. n->curpoint[[2]]],
  634.                        openlist = Append[ Drop[openlist, -1],
  635.                               { lastpoint[[1]],
  636.                                 lastpoint[[2]],
  637.                                 curpoint[[2]] } ];
  638.                        curpoint = {} ] ];
  639.  
  640.                  If [ Length[simplelist] > 0 &&
  641.                     Length[curpoint] > 0,
  642.                   nextpoint = First[simplelist];
  643.                   If [ SameQ[curpoint[[2]], nextpoint[[2]]-1] &&
  644.                        SameQ[curpoint[[1]] /. n->curpoint[[2]],
  645.                           nextpoint[[1]] /. n->curpoint[[2]]],
  646.                        AppendTo[ openlist,
  647.                          { nextpoint[[1]],
  648.                            nextpoint[[2]] - 1,
  649.                            nextpoint[[3]] } ];
  650.                        simplelist = Rest[simplelist];
  651.                        curpoint = {} ] ];
  652.  
  653.                  If [ Length[curpoint] > 0,
  654.                   AppendTo[openlist, curpoint] ] ] ];
  655.  
  656.         lastpoint = Last[openlist];
  657.         While [ SameQ[First[lastpoint] /. n -> lastpoint[[3]], 0] &&
  658.             ! SameQ[lastpoint[[3]], Infinity],
  659.             openlist = If [ SameQ[lastpoint[[2]], lastpoint[[3]]],
  660.                     Drop[openlist, -1],
  661.                         Append[ Drop[openlist, -1],
  662.                         { lastpoint[[1]],
  663.                           lastpoint[[2]],
  664.                           lastpoint[[3]] - 1 } ] ];
  665.             lastpoint = Last[ openlist ] ];
  666.  
  667.         openlist ]
  668.  
  669. (*  SetConvolutionDomain  *)
  670. SetConvolutionDomain[domain_] :=
  671.     Block [    {},
  672.         Unprotect[ $ConvolutionDomain ];
  673.         $ConvolutionDomain = domain;
  674.         Protect[ $ConvolutionDomain ];
  675.         domain ] /;
  676.     SameQ[domain, Discrete] || SameQ[domain, Continuous]
  677.  
  678. SetConvolutionDomain[bogusdomain_] := Message[SetConvolutionDomain::domain]
  679.  
  680.  
  681. (*  E N D     P A C K A G E  *) 
  682.  
  683. End[]
  684. EndPackage[]
  685.  
  686. If [ TrueQ[ $VersionNumber >= 2.0 ],
  687.      On[ General::spell1 ];
  688.      On[ General::spell ] ]
  689.  
  690.  
  691. (*  A L I A S E S  *)
  692.  
  693. PiecewisePlot = PlotList
  694. PiecewisePlot::usage = PlotList::usage
  695.  
  696.  
  697. (*  H E L P     I N F O R M A T I O N  *)
  698.  
  699. Block [    {newfuns},
  700.      newfuns =
  701.      { AutoCorrelation,        ConvertFromList,    ConvertToList,
  702.        CTPiecewiseConvolution,    DTPiecewiseConvolution, IntervalQ,
  703.        PiecewiseConvolution,    PiecewisePlot,        PlotList,
  704.        SimplifyList,        ValidIntervalQ };
  705.      Combine[ SPfunctions, newfuns ];
  706.      Apply[ Protect, newfuns ] ]
  707.  
  708.  
  709. (*  E N D I N G     M E S S A G E  *)
  710.  
  711. Print["Piecewise convolution rules have been loaded."]
  712. Null
  713.